function edgemap = get_eyelid(iris_location,leftpart,theta_range,isTop,debug)
    cx = iris_location(1);
    cy = iris_location(2);
    radius = iris_location(3);
    edgemap = edge(leftpart,'canny',[],2);
    edgemap(end-10:end,:) = 0;
    edgemap(1: cy - radius - 5,:) = 0;
    edgemap(cy + radius + 5:end,:) = 0;
    %[edgelist, ~] = edgelink(edgemap, 20);
    
    % Deal with the left part
    edgemap_upper = edgemap;
    
    [H,T,R] = hough(edgemap_upper,'RhoResolution',5,'Theta',theta_range);
    P  = houghpeaks(H,5);
    lines = houghlines(edgemap_upper,T,R,P,'FillGap',5,'MinLength',20);
    
    if isempty(lines)
        lines = houghlines(edgemap_upper,T,R,P,'FillGap',5,'MinLength',10);
    end
    
    if isempty(lines)
        lines = houghlines(edgemap_upper,T,R,P,'FillGap',5,'MinLength',5);
    end
    contrast = zeros(1,length(lines));
    
    if debug
%         figure(6); imshow(edgemap_upper); hold on;
%         for line_index=1:length(lines)
%             xy = [lines(line_index).point1; lines(line_index).point2];
%             %plot(xy(:,1),xy(:,2),'r');
%         end
    end
    
    [y_edge,x_edge] = find(edgemap_upper>0);
    
    for line_index=1:length(lines)
        xy = [lines(line_index).point1; lines(line_index).point2];
        minx = min(xy(:,1));
        maxx = max(xy(:,1));
        %figure(6); imshow(edgemap_upper); hold on; plot(xy(:,1),xy(:,2),'r');

        var_map_upper = zeros(size(edgemap));
        for j=1:length(x_edge)
            xi = x_edge(j); yi = y_edge(j);
            x2x1 = lines(line_index).point2 - lines(line_index).point1;
            x1x0 = lines(line_index).point1 - [xi yi];
            d = abs(det([x2x1' x1x0']))/sqrt(sum(x2x1.^2));
            var_map_upper(yi,xi) = d;
        end

        upper_left = edgemap_upper;
        var_map_upper = var_map_upper/max(var_map_upper(:));
        upper_left(var_map_upper > .05) = 0;

        upper_left(:,1:minx-5) = 0; upper_left(:,maxx+5:end) = 0;

        % Rotate the image
        if lines(line_index).theta > 0
            rotated_image = imrotate(leftpart,lines(line_index).theta-90);
            rotated_edgemap = imrotate(upper_left,lines(line_index).theta-90);
        else
            rotated_image = imrotate(leftpart,lines(line_index).theta+90);
            rotated_edgemap = imrotate(upper_left,lines(line_index).theta+90);
        end
            
        [y,x] = find(rotated_edgemap > 0);
        minx = min(x); maxx = max(x);
        miny = min(y)-10; maxy = max(y)+10;
        height = maxy-miny;
        bb = [minx miny maxx-minx maxy-miny];

        if debug
            h=figure(20);
            imshow(rotated_image); hold on; rectangle('Position',bb);
            plot([minx maxx],[miny+floor(height/2) miny+floor(height/2)],'r');
            pause(.1);
        end
        
        top = rotated_image(miny:miny+floor(height/2),minx:maxx);
        bottom = rotated_image(miny+floor(height/2)+1:maxy,minx:maxx);
        contrast(line_index) = exp(mean(double(top(:)))) - exp(mean(double(bottom(:))));
    end
    
    % Get the best line
    if isTop
        [~,mine] = min(contrast);
    else
        [~,mine] = max(contrast);
    end
    xy = [lines(mine).point1; lines(mine).point2];
    minx = min(xy(:,1));
    maxx = max(xy(:,1));
    if debug
        %figure(6); imshow(edgemap_upper); hold on; plot(xy(:,1),xy(:,2),'r');
    end
    result = zeros(size(edgemap));
    
    for j=1:length(x_edge)
        xi = x_edge(j); yi = y_edge(j);
        x2x1 = lines(mine).point2 - lines(mine).point1;
        x1x0 = lines(mine).point1 - [xi yi];
        d = abs(det([x2x1' x1x0']))/sqrt(sum(x2x1.^2));
        result(yi,xi) = d;
    end
    
    result = result/max(result(:));
    edgemap(result> .05) = 0;
    edgemap(:,1:minx) = 0; edgemap(:,maxx:end) = 0;